
rm(list=ls(all=TRUE))
sink("log_replication_WindowsRandom.txt") 




library(haven)
library(plyr)
library(sp)
library(foreach)
library(doMPI)
# create and register a doMPI cluster if necessary
cl <- startMPIcluster()
registerDoMPI(cl)

setwd("/scratch/user/fhollenbach/clusterRandomWindows/")

llCRS<-CRS("+proj=longlat +ellps=WGS84")


## let's just load the replication data
load("africa_shape_grid_windows.RData")
load("W_windows.Rdata")

########
conflict = read.csv("UCDP GED v.1.0-2011.csv")  



conf2008<-subset(conflict, Year==2008)



windowsStart <- round(seq(1,round(dim(conf2008)[1]/2),length.out = 45))
windowsEnd <- round(seq(round(dim(conf2008)[1]/2),dim(conf2008)[1],length.out = 45))


africa_shape_grid@data$spat_lag = as.numeric(mat%*% africa_shape_grid@data$conf2007_dum)

randomWindows = function(seeds,iteration,conf_data, shapefile_data, window, start, end){
  set.seed(seeds)
  conf_data$random.order = runif(dim(conf_data)[1],1,dim(conf_data)[1])
  conf_data <- conf_data[order(conf_data$High_est,conf_data$random.order),]
  
  sub = conf_data[c(start[window]:end[window]),]
  coordinates(sub)=~Lon+Lat
  proj4string(sub)=llCRS
  EventsInWindow = mean(sub$High_est)
  sub <-SpatialPoints(sub,proj4string=llCRS)
  points_sub2007 <-over(shapefile_data,sub,returnList=TRUE)
  for (j in 1:10674){
    shapefile_data@data$conf_window[j] <- ifelse(length(points_sub2007[[j]])>0,1,0)
    }
  mod1 = glm(conf_window ~ pre2000_count+bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + spat_lag,data=shapefile_data@data, family = binomial(link= "logit"))
  coefficients = rnorm(300,summary(mod1)$coef[[9,1]], (summary(mod1)$coef[[9,2]])^2)
  res = data.frame(rep(window,300), rep(iteration, 300), rep(EventsInWindow, 300), coefficients)
  return(res)
}

randomWindows_Best = function(seeds,iteration,conf_data, shapefile_data, window, start, end){
  set.seed(seeds)
  conf_data$random.order = runif(dim(conf_data)[1],1,dim(conf_data)[1])
  conf_data <- conf_data[order(conf_data$Best_est,conf_data$random.order),]
  
  sub = conf_data[c(start[window]:end[window]),]
  coordinates(sub)=~Lon+Lat
  proj4string(sub)=llCRS
  EventsInWindow = mean(sub$Best_est)
  sub <-SpatialPoints(sub,proj4string=llCRS)
  points_sub2007 <-over(shapefile_data,sub,returnList=TRUE)
  for (j in 1:10674){
    shapefile_data@data$conf_window[j] <- ifelse(length(points_sub2007[[j]])>0,1,0)
  }
  mod1 = glm(conf_window ~ pre2000_count+bdist1+capdist+pop2005+mnt+irri+gcppc00+cell_dum_07 + spat_lag,data=shapefile_data@data, family = binomial(link= "logit"))
  coefficients = rnorm(300,summary(mod1)$coef[[9,1]], (summary(mod1)$coef[[9,2]])^2)
  res = data.frame(rep(window,300), rep(iteration, 300), rep(EventsInWindow, 300), coefficients)
  return(res)
}



iterations = seq(1,250,1)
windows = seq(1:45)
toDo = expand.grid(iterations, windows)
names(toDo) = c("iterations", "windows")
toDo$seed = seq(1,dim(toDo)[1],1)
head(toDo)

dim(toDo)
random_high_results = foreach(i = 1:dim(toDo)[1], 
        .export = c("africa_shape_grid","conf2008","windowsStart", "windowsEnd", "toDo"),
        .combine = rbind,
        .inorder = FALSE
                              )%dopar%{
              result =  randomWindows(seeds = toDo$seed[i], iteration = toDo$iterations[i], conf2008, africa_shape_grid, window = toDo$windows[i], windowsStart, windowsEnd) 
              return(result)
             }
save(random_high_results, file = "random_high_results.rda")


random_best_results = foreach(i = 1:dim(toDo)[1], 
                              .export = c("africa_shape_grid","conf2008","windowsStart", "windowsEnd", "toDo"),
                              .combine = rbind,
                              .inorder = FALSE
                              )%dopar%{
    result =  randomWindows_Best(seeds = toDo$seed[i], iteration = toDo$iterations[i], conf2008, africa_shape_grid, window = toDo$windows[i], windowsStart, windowsEnd) 
    return(result)
  }
save(random_best_results, file = "random_best_results.rda")



closeCluster(cl)
mpi.quit()

sink()

